Response and fluctuations of a two-state signaling module with feedback 
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We study the stochastic kinetics of a signaling module consisting of a two-state stochastic point 
process with negative feedback. In the active state, a product is synthesized which increases the 
active-to-inactive transition rate of the process. We analyze this simple auto-regulatory module 
using a path-integral technique based on the temporal statistics of state-flips of the process. We 
develop a systematic framework to calculate averages, auto-correlations, and response functions by 
treating the feedback as a weak perturbation. Explicit analytical results are obtained to first order 
in the feedback strength. Monte Carlo simulations are performed to test the analytical results in the 
weak feedback limit and to investigate the strong feedback regime. We conclude by relating some 
of our results to experimental observations in the olfactory and visual sensory systems. 
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I. INTRODUCTION 

Signal transduction in biological cells refers to the set 
of processes by which the cells receive and process infor- 
mation from the environment. From the realization of 
the extra-cellular stimulus to its final recognition by the 
cell/organism, there is, in general, a complex biochemi- 
cal reaction network, commonly referred to as the signal 
transduction pathway. Faithful transmission of informa- 
tion through the pathway may require amplification and 
adaptation, which is often accomplished through posi- 
tive (amplification) and negative (adaptation) feedback 
mechanisms built into the reaction network. 

External and intrinsic noise can potentially limit the 
faithful transduction of a 'signal' (information encoded 
in the intensity of the external stimulus and its spatio- 
temporal variation) [l], Since cellular biochemical 
networks often function with a limited number of react- 
ing molecules, noise is always present, and an important 
question is how the cell physiology maintains its robust- 
ness in spite of the randomness in the underlying molec- 
ular events Q. In the context of signal transduction, ex- 
ternal noise refers to the randomness in the signal (fluc- 
tuations in the local concentration of the extracellular 
ligand etc.), whereas intrinsic noise is a collective effect 
produced by the inherent stochasticity of the transduc- 
tion mechanisms itself, such as random opening and clos- 
ing of the ion channels, fluctuations in reaction rates and 
so on 0, i, i, 0, 1, la [Mini- Although randomness 
and noise are beneficial in some instances (population 
heterogeneity being a well-known example [12j|). by and 
large, cells have evolved the means to control biochemical 
noise through regulatory mechanisms (such as negative 
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feedback) to ensure the robustness of biochemical net- 
works [13]. Indeed, in the context of electrical circuits, 
negative feedback is a well-known noise-reduction mech- 
anism used in devices such as amplifiers and oscillators 
(see [13] for a review and important references). 

Complex signaling pathways can often be pro ductively 
viewed as consisting of recurring modules jlfj . flfl ]. A 
module is typically made up of multiple species of inter- 
acting molecules acting together with a specific function. 
A signal transduction pathway may thus be described 
in terms of a series of modules, interacting with each 
other 17]. In very general terms, a module receives a 
signal (from the extracellular environment, or from an- 
other module) and transmits it, possibly in a different 
form (e.g. chemical signals converted to electrical pulses). 
Each module can function to a certain extent in isolation, 
but interaction among modules and their coordination 
are crucial for carrying out all the complex functions of 
the cell. 

In this paper, we analyze a signaling module based on 
a single protein that switches between active (A*) and in- 
active (A) states. In the active state, a certain molecular 
species C is produced with a fixed rate. The accumula- 
tion of C increases the A*-^A transition rate (Fig. [l}, 
leading to negative feedback on the production of C. In 
addition, C is removed at a fixed rate independent of the 
activation state. As an example, Fig. [2] illustrates the 
modularity of signal transduction and the occurrence of 
the module defined in Fig. [TJ for the specific case of ver- 
tebrate olfactory sensory neurons (reviewed in fill [l9l]). 
(A closely related pathway also operates in cone photore- 
ceptors in the retina.) 

In the work presented, we develop a systematic analyti- 
cal framework to study the dynamics of the simple auto- 
regulatory module shown in Fig. [TJ Without the feed- 
back, the activation and deactivation kinetics would be 
described by a two-state Poisson process. In the presence 
of feedback, the deactivation rate, at any given time, be- 
comes dependent on the history of transitions, and makes 
the effective two-state kinetics non-Markovian. We use 
a path-integral technique to evaluate the effects of feed- 
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FIG. 1: The simple two-component auto-regulatory signaling 
module with negative feedback studied in this paper. The 
product C (with concentration c), generated with rate J in 
the active state A*, enhances the deactivation rate R-. C 
is removed with a rate A independent of the activation state 
of A. An input signal encoded in a temporal change R+(t) is 
transduced into an output signal c(t). 

odorant 



FIG. 2: A schematic representation of the modularity of sig- 
nal transduction in olfactory sensory neurons. (The arrows 
with + sign indicate enhancement of the reaction rates, arrows 
without signs indicate the direction of the reactions in the con- 
ventional sense.) After binding of an odorant (external stim- 
ulus) to the receptor R, the enzyme adenylate cyclase (AC) is 
activated and produces the 'second messenger' cAMP, leading 
in turn to the opening of cyclic-nucleotide-gated (CNG) ion 
channels and to the depolarization of the cell membrane. Two 
negative feedback mechanisms (reviewed in 20]) are shown 
in the figure: (i) cAMP-activated kinase phosphorylates (and 
therefore deactivates) the receptor R, (ii) Ca 2+ ions bound 
to calmodulin inhibit the Ca 2+ channels, which is crucial for 
adaptation to repeated odorant stimuli [2ll . I22I ]. The signal 
transduction pathway can be viewed as consisting of two mod- 
ules of the type shown in Fig. [T] connected in series. 



back perturbatively, which enables us to derive the noise 
from the dynamics of the system. This is in contrast 
to stochastic analysis of biochemical networks based on 
Langevin equations [23| with additive white noise (see, 
e.g. MM)- 

We obtain explicit analytical expressions for the differ- 
ent statistical averages as well as correlation and response 
functions of the module in the limit of weak feedback. We 
also perform Monte Carlo simulations to test the analyt- 
ical results in the weak feedback limit and to investigate 



the strong feedback regime. The path-integral formalism 
we develop is similar to the one described in [25|, and 
the system considered is the same as in Q. The calcu- 
lations in however, focus on steady state probability 
distributions rather than temporal correlations. 

The paper is arranged into four sections. Sec. [IT] pro- 
vides a brief description of our simple model of a stochas- 
tic signaling module with feedback. All important quan- 
tities are introduced and the notation is defined. In 
Sec. IIIIl we present and develop the path-integral for- 
malism. In Sec. IIV1 we compute the different Green's 
functions using this formalism to first order in feedback 
strength. Sec.[V]is devoted to the calculation of the aver- 
ages and correlation functions, and Sec. I VII contains the 
calculation of the response functions. Sec. IVIll presents 
the results of Monte Carlo simulations and their compar- 
ison to analytical results. Our conclusions and outlook 
for further extensions are presented in Sec. IVIIII Four 
appendices detail some of the analytical calculations. 



II. TWO-STATE MODEL: GENERAL SET-UP 

We describe our model of a two-state signaling mod- 
ule using the example of an ion channel connected to a 
small cellular compartment (Fig. [3]). The channel on the 




FIG. 3: A channel in the open state permits the entry of ions 
(rate J) into a small compartment of volume V. The ions, 
once inside, increase the closing rate of the channel, and are 
also removed from the compartment at a constant rate A. 

cell membrane exists in one of two conformational states: 
open, when the ions enter through the channel, and closed 
when this is not allowed. The state of the channel at time 
t is represented by S(t), with S = 1 being the open state 
and S = 0, the closed state. The — > 1 transition takes 
place with a rate and the reverse transition 1 — > 
occurs with a rate Let us denote by c(t), the ion 

(which in the case of the olfactory signaling pathway is 
Ca 2+ ) concentration inside the compartment at time t. 
The kinetics of c is described by the equation 



dc 



S(t)-Xc(t). 



(1) 



In the above equation, J represents the molar current 
of ions entering the cell through the channel in the open 
state, A represents the total rate of removal of ions by 
membrane pumps, and V is the volume of the compart- 
ment. The channel kinetics is specified by the opening 
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rate R+ and closing rate We assume the existence 
of a negative feedback corresponding to a rate i?_ that 
increases with increasing ion concentration. If the effect 
of the feedback is weak, one may expand to linear order 
and write 

R-{c)k R°_+ac(t), (2) 

where R°_ is the closing rate of the channel when no ions 
are present in the compartment and the coupling param- 
eter a specifies the feedback strength. We further as- 
sume that, in general, an external signal is received by 
the module as a change in the opening rate (due, e.g., to 
the binding of extracellular ligands to the channel). We 
may write 

R + $)=R° + +$(i), (3) 

where RP, is the intrinsic opening rate of the channel in 
the absence of external stimulus, and 4>(t) defines the 
stimulus. Eqs. [TJ [5] and [3] specify the dynamics of the 
problem. 

It is now convenient to adopt a dimensionless formula- 
tion of the problem. The inverse of the intrinsic closing 
rate R°_ is chosen as the unit of time and the ratio J /(XV) 
(i.e. the maximum achievable ion concentration) is the 
unit of c. Using these two quantities, we scale all the 
other parameters, and the complete list of dimensionless 
parameters is given below: 

XV aJ/V x A (j> 

c = —~ c; a = Tk^ ; X = W ; * = W' 

r+ = |±; r_ = -^ = l + ac; t = R?_i. (4) 

The dynamical equation for c(t) is then simplified to 4| = 
X(S(t) — c(t)), with the solution 

c(t) = A f e-^-^S^dt 1 . (5) 

J — oo 

Fig. H] shows two typical time evolution curves of S(t) 
and c(t). Note, from Eq. [5l that in the limit A — * oo, 
c(i) ps S(t), which is illustrated in Fig. 0J In this limit, 
the Ca 2+ is drained from the compartment as soon as 
the channel closes, and consequently, the time evolution 
of c(t) closely follows the channel state. In this sense, A is 
an important control parameter in our model, and deter- 
mines how much the dynamic characteristics of c(t) (in- 
cluding fluctuations and response functions) are tied to 
the corresponding quantities for the channel state. These 
aspects will be discussed more in Sec. IVland later. 

The combined set of Eqs. [T][3] describe the kinetics of 
S(t), which may be characterized using, for instance, the 
n-point functions (n — 1,2,...): (S(t Q )...S(t n -i)). The 
corresponding functions for c(t) are then easily computed 
from these using Eq. [5] 

(c(* )...c(*„-i)) = X n e- x ^+- +t ^ x 
T K- r"<_ 1 e A ^ + - +t «-^<5(t' )...5(4_ 1 )).(6) 

J — oc J —oo 
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FIG. 4: The figure illustrates the kinetics of channel openings 
and Ca 2+ concentration for two typical runs, with r+ — 0.5 
and A = 5 (top and middle) and A = 0.5 (top and bottom) in 
the absence of feedback (a = 0). Note that for larger A, c(t) 
almost follows the step-like kinetics of S(t). 



The angular brackets (• • • ) represent statistical averaging 
over the different temporal histories of the process. In 
particular, when the external signal is time-independent, 
the system reaches a steady state, where the one-point 
functions (S) and (c) are constants, while the auto- 
correlation and cross-correlation functions 

C 5 (t)= lim (S(t )S(t + t)) -(S) 2 , 

tQ — »oo 

C c (t) = lim (c(t )c(t + t)) - (c) 2 , 

to — >oo 

C Sc (t) - lim (S(t )c(t + t)) - (S)(c), (7) 

tQ — >OG 

are stationary in time. In particular it also follows from 
Eq. [5] that 

(c) = (S) (8) 

in the steady state. We may also define the power spec- 
tral densities for the fluctuations in S and c from the 
stationary auto-correlation functions, for example, 

f°° X 2 
P c (uj) = 2 / C c (t) cos(ujt)dt = 2 Ps(u), (9) 

Jo A + UJ 

where the relation between P c and Ps follows from Eq. [S] 
Note that for to <C A, P c ~ Ps, whereas when uo ^> A, 
P c ~ Lu~~ 2 Pg. This has an interesting physical interpre- 
tation: when A is small over the time scales of interest, 
Ca 2+ kinetics is slow, and this effectively suppresses the 
power-spectrum at larger frequencies compared to the 
case of large A. 

The response of the system to a time-dependent exter- 
nal perturbation 4>(t) is {S{t)) <t > - (S) and (c(i))* - (c), 
where the superscript tf> indicates that the function needs 
to be evaluated in the presence of the perturbation. 
When the external signal is weak, i.e. <f>(i) <C r+, the 
response is characterized using the linear response func- 
tions Xs{t) and Xc(t) defined through the following rela- 
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tions: 



(10) 



Using Eq. [5j the following relation between the response 
functions is also easily proved: 

Xc(t) = X f xs(t')e- x ^dt'. (11) 
Jo 

In the following sections, we present explicit calcula- 
tions of the correlation and linear response functions of 
the model in the presence of feedback. Most of the spe- 
cific results in this paper are derived for a constant r + . 



III. PATH-INTEGRAL FORMULATION 

We start by introducing the propagator Uy (io, Co; t, c) 
with i,j = {0,1} that gives the probability density to 
find S(t) — j and c(t) — c given that S(to) — i and 
c(to) = cq. The mean open fraction of the channel in 
steady state can be written as 

(S) = f n a (-oo,co;0,c)(ic (12) 
Jo 

and is independent of i and Co- We define two different 
mean values of the concentration c, namely the mean con- 
centrations during periods when the channel is closed and 
open, respectively, starting from a special initial state 
S(t = 0) = 0, c(t = 0) = 0: 



(c(t)) = / n i0 (0,0; t,c)cdc, 
Jo 

(c(*))i = / H<i(0,0; t,c)cdc 



(13) 



The auto-correlation function of S in steady state is 

(S(0)S(t))= I del dciIIi 1 (-c»,e»);0 ) c 1 )nii(0,ci;t,c). 
Jo Jo 

(14) 

The linear response function can be written using a step 
function r + — > r + + 0(f) (with = <f>o&(t)) as a stim- 
ulus: 



E 



_9 _a_ 

Si <9<^>o 

n#(0,ci;t,c) 



=o , 



i ,1 

dc / dcillij(— oo, Co; 0, ci)x 
'o 



(15) 



In order to formulate a path-integral representa- 
tion of the propagator, we define the functional 
r Pij[to,t;{T k }%L 1 ;c ;N] of S(t) with i > t . It repre- 
sents the probability that a certain time evolution of the 



system from S(to) = i and c(to) — c to S(t) = j is re- 
alized. This time evolution is characterized by the set of 
times {rk}^ =1 at which the state switches from S = to 
1 or from S = 1 to 0. The number N represents the total 
number of such state changes in the time interval [to,t] 
and is even if i = j and odd otherwise. The propagator 
Hy (to, Co; t, c) can be expressed as a sum over all possible 
realizations of S(t) as: 

IIy(to,co;t,c) = 53 / 'W i j[h,t]{n}k=i><%'> N \S[c(t)-c], 

N J 

(16) 

where fVr = J* t * dn dr 2 ... f* N i dr N . 



It is convenient to introduce the reduced propagator 
or Green's function 



Gij(to,co]t) = / Ilij(to,co]t,c)dc, (17) 
Jo 

which can be written as 

G ij (t ,c ;t)=J2 I ©rPy[to,t;{T fc }f =1 ;c ;Ar]. (18) 

JV 

By definition the following relations hold: 

Gqx + Goo = 1 = Gio + Gii. (19) 

Note that in the absence of feedback, the func- 
tional Voi and V\\ do not depend on c, and the 
history-dependence is thus absent. In this case, the 
expression for the n-point functions reduces to the 
product of Green's functions: (S(to)...S(t n -i)) = 

G^ ) (-oo,i )G^ ) (t ,ii) • • •G^° 1 ) (t n _ 2 ,i„-i), where the 
superscript (0) denotes the absence of feedback (valid for 
the rest of the paper and all quantities) [iij . As special 
cases, the one- and two-point functions in the absence of 
feedback are given by 

(S(to)) = (S) = G °i ) (-oo,i ) ) 



(S(t )S(t)) = G^(-oo,t )G^(t ,t), (20) 
when a = 0. 

We now write down the functionals Vij in the explicit 
form. The expressions are presented only for Vqq, exten- 
sion to the other cases are straightforward. Furthermore, 
since the different Green's functions are related (Eq. [P9l 
and Eq. [37] later), it is enough to compute one of them. 
It is instructive to start with a — 0, where the transi- 
tions — > 1 and l—>0 are characterized by the (time- 



Co), 



invariant) rates r + and r_ 



1, respectively (in 



terms of scaled time). Let us consider a certain history 
of the process such that the channel is closed at t = 
and t, but changes state an even number N times in be- 
tween, at times n, ra, tjv- First consider the interval 
[0, 7~i]. The probability that the channel remains in state 
until time T\ has the form Po(0,Ti) = e~ r+Tl . Simi- 
larly, the probability of the channel to remain in state 1 
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during the time interval [ti,t 2 ] is -Pi(ti,t 2 ) = e ' T2 Tl ^. 
The (differential) probability for the whole process is 



Ho [o, *; {n}g 1 ;N]dT 1 dT 2 ...dT N = 



of the dimensionless feedback strength a as follows: 



Poo[0, *;{«*},{*$};<*; 2m] = ripe 

2m 2m— 1 



YldnPo(0,T 1 )P 1 (T U T 2 )P (T2,n)---P (T N ,t). (21) \z^< T i)~l^ / drc ( r ) 

7=1 \<=a i=i ^ 



1 + a x 



+ 0(a 2 ). (26) 



We now substitute for pj and Pi , and make the trans- 
formation from flip-times {t{\^ =1 to time-interval vari- 
ables {tj,t'j}™ =1 with m = AT/2, where the {£.,} denote 
the times, when the channel is closed and the {tj} when 
it is open (see Fig. [5]). The relation between the flip-times 
and the time intervals is summarized below: 



ti = h, 

1/2 i/2 

Ti = 2_, tl + t'i even i > 2, 

i=i i=i 

(j + l)/2 (i-i)/2 

t j = ti+ J2 *J odd j > L 



(22) 



S 

1 + 







~ n 



A 





















r 2 r 3 r 4 



T2m-1 T2 n 



FIG. 5: A schematic diagram of time evolution of S(t), show- 
ing the time flip variables n versus the interval variables ti 
and t\. 

The Jacobian for the transformation is J = 1, as fol- 
lows from Eq. [22] The probability functional thus has 
the form 



P$[0, t;{ti}, {4}; 2m] 



r m e -f O o(0,t;{t-};m)^ ^3) 



with the weight factor Poo given by 

m 

P oo (0, t; {Q; m) = r + t + (1 - r+) ]T 4 (24) 



»=i 



In the presence of the feedback-coupling to c(i), the 
expression for Pqo is modified to 



Y[ {l + ac(n))e 



Voa[0,t; c ; 2m] =r™ x 

E„2m-1 
/" J + 1 drc(r) 

3 = 1 3 



i=2 



where the prime (double prime) on the product- or sum- 
symbol indicates that the running index is always even 
(odd) . The above expression may be expanded in powers 



Eq. [26] is the basis for more specific and detailed calcula- 
tions to follow in the next sections. 



IV. COMPUTATION OF GREEN'S FUNCTIONS 
A. Computation of Goo 



For instructive purposes, we first show how the Green's 
function Gqq (0, t) in the absence of feedback is computed 
using our formalism. In this case the exact answer is 
easily found by solving its rate equation 



ftG$=r--(r_+r + )G&\ 



(27) 



with T- = r°_ = 1 and the initial condition G °q (0, 0) = 1. 
The solution is 



G, 



(0,t) = —?—(l+r + e-( 1 + r ^). (28) 
1 + r+ V / 



This Green's function can also be computed using the 
path-integral technique outlined previously, more specif- 
ically from Eq. [18] using Eqs. [23] and [24] For simplicity, 
let us put to = 0, and define the Laplace transform of 
the Green's function G$(s) = J °° G$ (0, t)e~ st dt. The 
calculation is most easily done using the generalized con- 
volution theorem presented in Appendix A fEq. lAip . Us- 
ing this technique, the path integral in Eq. [TS] becomes a 
product in the s-space, and the result is 



G { o°o\s) =9(s + r+), 
l °° /lY 

with *«=-£(-) 



s + l-r. 



(29) 



After summing the geometric series in Eq. 1291 we find 

that G ° \s) = (s + I)/ [s(s + 1 + r+)], which, upon in- 
version, gives Eq. |2"51 

Let us now extend the previous calculation to include 
feedback, and compute Goo(0,co;i) to first order in a: 

G O o(0,c o ;i) = G$(0,t) +aG$(0, c ; t) + 0(a 2 ). (30) 

(25) From Eq. [55] and Eq. UM we find that G^ (0, c ; t) de- 



pends only linearly on Co: Gq ^(0, Co;t) = Gqq (0, Co 
0; t) + cof(t). Putting together the co-independent terms, 
we write 

G 00 (0, c ; t) = G oo (0, c = 0; t) + ac f(t) + 0(a 2 ) : (31) 
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where co/(i) is simply the 0(a) term in Goo(0, c ; t), 
when the channel evolves in a situation with J = so 
that r_(t) = 1 + acoe~ xt . Using this r_(t), the function 
/(i) is easily determined by solving the rate equation 
Eq. [271 with the initial condition G oo (0, c ;0) = 1. The 
result is 



/(*) = ^ 



'"' 1 1 c -(l+r+)t ! 1 e -(l+r++A)f 

A \ 1 + r + - A 1 + r+ ' 



A 



-At 



(l + r+)(l + r+-A) 



(32) 



The function Goo(0,co = 0;i) itself may be expressed 
in terms of three integrals as follows: 

(oo 
1 + l>?[/o(*;m) 
m=l 

+a(/ 1 (t;m)-I 2 (*;m))]J + 0(a 2 ), (33) 

where 
I (t;m) 



h(t;m) 
h(t;m) 



2m 
2m-l T 



(34) 



Note that we have introduced the compact nota- 
tion JVT = Jl dh Sf h dt[ Jo"' 1 "*' 1 dfe... JT*-"— *~ dt ' m 
for the integration measure. The explicit calcula- 
tions are done most conveniently in terms of Laplace- 
transformed variables, defined in the standard way: 
f(s) — J °° f(t)e~ st dt. The previous Eq. [33] becomes, 
in this notation, 



G o(c =0;s) = (s + r+) 1 

OO 

+ X>+ [lo(s + r+;m) 
+a(l 1 (s + r + -,m)-l2(s + r + ]m)) +0(a 2 ). (35) 



The calculation of the integrals Iq, I\ and I2 is done in 
Appendix El We omit further details, and present only 
the final result: 



G oo (0, c = 0; t) = G$ (0, t) + aG$ (0, c = 0; t) = g£ } (0, t) + a 



(l+r+Y 



r+ + A 
1 + r+ + A 



-A/ 



r+(l + >+) 
(l + r+-A) 2 



,-(X+r++A)i_ 



l + r + 
A(l + r+ + A) 



(1+r+)t f (l + r+)(A-l)t r 2 _-(A-l) 3 + r + (2-3A + 2A 2 ) 



1 + r 1 



A(l + r+ - A) 2 



0(a 2 ). (36) 



B. Relation between Gn and Goo 

We will now show that there is a non-trivial relation 
between Gn(0, en = 0; t) and Gno(0, en = 0; t), as follows: 

Gn(0,c = 0;t) =G 01 (Q,c = Q;t) 

^Gqo(0,c = 0;t) 
+ ftGoo(0,co = 0;t)|t=o' 1 j 

which will be shown to be true up to 0(a). To prove 
this, let us start with the case a = 0. Then, the following 



relation is true for any arbitrary < t' < t: 

G<°>(0,t) = G$(0,t')G$(t',t) + G^(0,t')G^(t',t). 

(38) 

Let us now take the limit of t' — ► 0, and use the 
Taylor expansions G$(0,t) = 1 + t ' aG °^ 0>t) \ t=0 + • • • , 

Gf? (f , t) = Gf { t ~ t>) = CI? (0. «) - + " " " as 

well as the condition G^ = 1 — Gqq . After substituting 
these back into Eq. [351 we arrive at Eq. [37] 

In the presence of feedback, Eq. [33Jis not true anymore 
because of the explicit history-dependence. However, us- 
ing expansions of the two Green's functions Goi(i',c';i) 
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and Gii(i', d; t) similar to Eg. [3T1 

G 01 (t', d; t) = G i(t', d = 0;t) + ad ' h{t - f) + 0(a 2 ), 

G n (t', d; t) = G n (t', d = 0;t) + ad f 2 (t - t') + 0(a 2 ), 

(39) 

Eq.[37]can be shown to be valid also for a ^ (to 0(a)). 
The proof as well as the functions fi(t) and f 2 (t) are 
given in Appendix IB1 



c ~ 1 in the open state of the channel. The mean open 
fraction, therefore, is given by 



(S) 



l + r++a l + r+ (l + r+) 



- + Q(a 2 ). (44) 



It is easily verified that this expression is the limit of the 
one given in Eq. [40] to 0(a) for A ^> r+ and A ^> 1. 

Our perturbative result for the mean fraction of open 
channels (Eq. 20]) is contained in the results of @ • 



V. AVERAGES AND CORRELATORS 



B. Auto-correlation functions 



A. Averages and comparison to mean-field analysis 

The mean fraction of open channels in the steady state 
is now easily found as (S) = lim^oo [1 — Goo(0, cq = 0;t)], 
and the result is 



(S) = (c) =- 



r+ 



A 



(l + r + )(l + r + + A) 



Q{a 2 



(40) 



For comparison, we may also compute the same quan- 
tity using the mean-field approach. In general, the steady 
state obeys the relation 



r+(l - (S)) = (r-S) = (S) + a(cS). 



(41) 



In the spirit of mean-field analysis, we now assume that 
the fluctuations in S and c are independent, (cS) — 
(c)(5). Making use of the equality (c) = (S) in Eq.HTI 
(valid in the steady state) , one obtains a quadratic equa- 
tion for (S). The solution is 



(S)mf = 



v / (l + r+) 2 +4ar+- (l + r+) 



2o 



(42) 



A small- a expansion of the square root gives 



(S)mf = 



l + r+ (l + r+Y 



+ 0(a 2 ). (43) 



The mean-field treatment can be expected to be valid 
if the calcium dynamics is slow compared to the channel 
dynamics, as in this case the calcium concentration ex- 
hibits only small fluctuations around its mean value and 
can be replaced by a constant (as will be seen explicitly 
in Eq. Iij2")) . Indeed, the full result in Eq. |4H1 reduces to 
the mean- field result in Eq. [43] when A«r + and A <C 1. 
For all other values of A, the mean-field analysis under- 
estimates the effect of feedback. 

A second limiting case of interest is the situation where 
the pump rate is large: A — > oo, in which case, from Eq.(5] 
we find that c(t) « S(t) at all times t. In this case the 
closing rate may be approximated as r_ w 1 + a since 



To compute the two-point function of S, we put the 
definition of the reduced propagator (Eq. [17]) into Eq. [T4l 

(S{0)S(t))=[ dc 1 n i i(-oo,c D ;0,ci)G! u (0,ci;t). (45) 



Using Eq. [39] together with Eqs. [P2l and [T3l we arrive at 

(S(O)S(t)) = (5)Gu(0,cd ^0;t)+a(c) 1 f 2 (t)+O(a 2 ), 

(46) 

where Gn(0, cq = 0;t) is given by Eq. [37] Goo(0, c = 
0;t) is given by Eq. [36] fyft) is given by Eq. IBll and (c)i 
is the steady state value of the restricted average (c(t))i 
(Eq. [T3"]) and calculated in Appendix [C] Using Eq. 0(5] 
and the steady state average in Eq. [40] we find the auto- 
correlation function of S (Eq. [7]) to be 



Cs(t) = 



,-(l+r+)t 



a (^B\e 



-(!+»•+)*. 



(l + r+) 2 

de' xt + D ie -( 1+r + +A )* + E 1 te' {1+r+)t ^ + 0(a 2 ).(47) 

The coefficients B\ , C\ , D\ , E\ are given in Appendix iDl 
We observe that the first order correction term intro- 
duces two new time scales. Furthermore, the feedback 
introduces a term non-monotonic in time, whose sign de- 
pends on the relative values of A and r+. For A < 1 
and A > 1 + r + , E\ > 0, whereas for intermediate values 
1 < A < 1 + r+, Ei < 0. This term is the first indication 
that the feedback term introduces qualitative differences 
in the decay of the auto-correlation function. 

The corresponding correlator for c(t) is computed from 
Eq. [6] using the two-point function in Eq. 1461 The result 
is 



C c (t) 



r + A 



(e- At (l + r+)- 



(1 + r + ) 2 [(l + r + ) 2 — A 2 ] 
Ae-( 1+r +>* ) + a (B 2e -( 1+r +) t + C 2 e~ xt + 
D2e -(i+r + +\)t + E2te ~(i+r + )t + F 2te -^ + 0( Q 2 ).(48) 

The coefficients B 2 , C 2 , D 2 , Eq, , F 2 are given in Ap- 
pendix [D] Similar to the previous case, the first order 
correction term has introduced two new time scales, and 
unlike the previous case, there are two non-monotonic 
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terms in time. From Eg. 1481 and using Eq. [9] we now find 
the power spectrum for c-fluctuations to O(a): 



_A 2 



2a 



Pc ^ 2 (l + r + )(^ 2 + A 2 )[co 2 + (l + r + ) 2 ] 
B 2 {l + r+) C 2 X D 2 (l + r + + X) 



Lu 2 + {l + r+) 2 



+ A 2 lu 2 + (1 + r+ + A) 2 



E 2 [(l 



r_i_ r - w 



[(1 



f 2 (A 2 -^ 2 ) 
(A 2 + oj 2 ) 2 



0(a 2 \A9) 



For large oj, P r .(uj) ~ w~ 4 , which is true with and with- 
out feedback [45|]. However, for small frequencies, i.e., 
lu <C A in particular, the power spectrum for c effec- 
tively decays as lo~ 2 because large A attenuates the high- 
frequency (but still less than A) noise. 

From Appendix [Dj we see that F 2 < always, whereas 
E 2 < for 1 < A < 1 + r+ and E 2 > for A < 1 
and A > 1 + r + . The interplay of the time scales A 
and 1 + r + gives rise to non-trivial time-dependence in 
the auto-correlation function. This results in a variety of 
possible time courses, examples of which can be observed 
in numerical simulations presented in Sec. IVII1 

From the correlators Eq. 27J and Eq. 05] the fluctua- 
tions are also determined easily. We define the mean- 
squared fluctuations through 



(5S) 2 =(S 2 )-(S) 2 = (S)(1-(S))=C S (0), 
(6c) 2 =(c 2 )-(c) 2 = C c (0)- 



(50) 



The channel state fluctuation has the following char- 
acteristic: Since (SS) 2 = (S)(l — (S)), its maximum is 
always j, which is reached when (S) = ^, independent 
of feedback. The RMS fluctuations, to first order in a 
are given by 



l + r+\ 2 {1 + r+)(l + r+ + X) J 

c / r +A ( \ a 

dC ^(l + r + ) 2 (l + r + + X) [ 1+ 2 X 

r 3 + + r% (3A - 2) + r+(2A 2 - 4A - 3) - A(2A + 3) 
(l + r+)(l+r + + A)(l+r+ + 2A) 
+ 0{a 2 ). (51) 



For the channel state, we see that the first order correc- 
tion term simply reverses sign at the point r + = 1 (when 
the intrinsic opening and closing rates are the same) . The 
fluctuations are enhanced by feedback when r+ < 1 and 
suppressed when r + > 1. 

The calcium fluctuation can likewise be increased or 
decreased by feedback, depending on the values of A and 
r + . The precision of output of the signaling module is 
given by the relative fluctuation (normalized standard 



deviation) in Ca 2+ . Using Eqs. l40l and [511 we obtain 



8c 
Jc) 



X 



(r++A)(- 



A) 



2A) 



2(1 



2A)(1 



A) 



+ 0{a 2 ). (52) 



Clearly for A > \ the relative fluctuation is increased 
for any value of r+, while for sufficiently small A it is 
increased when r + is large and decreased when r+ < 1. 
The precision of output can thus be improved by feedback 
when the input is weak (r + is small). 



VI. RESPONSE FUNCTIONS 

A. Computation of xs and \c 

When a step-stimulus r+ + (f>o9(t) is applied to the 
system, the response function (Eq. 1 1 5[) in S can be writ- 
ten using the definition of the Green's function (Eq. [T7[) 
together with Eq. [351 as 



, . d d 
Xs{t) = didf 



(S)G?i 5 (0,0;i) + (l-(S))G^(0,0;i) 
+a (ff°(t)(c) + ft(t)(c)i)] + 0(a 2 ), (53) 



where 4>o appearing as superscript indicates that r + is 
to be replaced by r + + (j) . Putting in the two Green's 
functions (Eqs.l36land l37|) together with Eqs. IB II and the 
restricted calcium averages (computed in Appendix [C|l. 
the final result is 



xs(t > 0) 



1 + r 4 



C 3 e~ xt + J D 3 e- (1+r + +A)i + E 3 te-^ 1+r ^ t j + 0(a 2 )(54) 

where B 3 ,C 3 ,D 3 ,E 3 are given in Appendix ID] The co- 
efficient i?3 of the non-monotonic term is negative for 
A < 1 and A > 1 + r+ , and positive for 1 < A < 1 + r + . 
The contribution from these terms makes the decay of 
the response function faster (relative to a = 0) for small 
and large A and slower for intermediate A. 

The response function Xc(t) is now calculated from 
Eq. [53] using Eq. [TT] yielding to first order in feedback 



Xc(t > 0) 



A 



( 



e -Ai _ e -(l+r+)t 



(l + r+)(l + r+-A) 



+E 4 te~ {1+r + )t +F 4 te 



-At 



+ 0(a 2 ), 



(55) 



where F 4 < 0, while E 4 < for A < 1 and A > 1 + r + 
and E 4 > for 1 < A < 1 + r + . 
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B. Violation of the fluctuation-dissipation theorem 

Having computed both the correlation and response 
functions for S(t) and c(t) separately, we will now show 
that these quantities violate the fluctuation-dissipation 
theorem (FDT) which holds for systems in thermal equi- 
librium. In Fourier-space, the FDT reads (with the 
Boltzmann constant ks and the temperature T) [26j 



u>C(w) — 2fcsTImx(w), 



(56) 



if the correlation function C{ui) and the response function 
x(w) refer to conjugate stimulus and response variables. 

To relate to the system considered in this article, an 
energy scale AU is introduced by writing the ratio of 
the transition rates (with c = 0) as R+/R°_ = r + = 
exp(— /3A.U), where /3 is the inverse temperature and AU 
is the energy difference between open and closed channel 
states. Let us now introduce a 'field' h which changes 
the energy gap (analogous to an external magnetic field 
in the context of magnetic models): AU — > AU + h, 
so that the dimensionless rate r+ is changed to r' + = 
r + e~P h sa r+(l — (3h) (where we have assumed that h is 
small). The relation between the energy change h and 
the (dimensionless) change in the flipping rate <p (which 
we defined earlier as the stimulus (Eq. [3])) is therefore, 
<j> = —f3r + h. Using this to transform our (dimensionless) 
response function xs(t), Eq. EH can be written dimen- 
sionless as 



uC s (u) = -2r + ImxsM, 

and analogously for c. 

In the absence of feedback, we have 



(57) 



2r+ 



(l + r + )[(l + r + ) 2 +c^ 2 ]' 

XsM = 77- TfA (58) 

(1 + r+)(l + r+ + iLu) 

from which we find (unsurprisingly) that the FDT holds 
for S(t). However, for c(t) the analogous results are (with 
a = 0) 



C c {lo) 



2A 2 r+ 



(l + r + )[(l + r+) 2 +w 2 ](A 2 +w 2 )' 

* cM = (l + r + )(A + iw)(l + r + +iu,)' (59) 

implying that FDT is violated for c(t), i.e. the dynamics 
of c cannot be represented by an equilibrium system with 
the same temperature T that governs the dynamics of S. 

In order to characterize this violation quantitatively, it 
is convenient to introduce an 'effective temperature' T cS , 
defined as the ratio between the left and right hand sides 
of Eq. US multiplied by T 0: 



2r + Imxs(u) 



(analogously for c). T is the 'real' temperature as intro- 
duced above. For a = 0, from Eq. [58] and Eq. [59j we 
then find that Tf/T = 1 for S (i.e. FDT is fulfilled), 
while Tf/T = A/(l + A + r+) < 1 for calcium (i.e. FDT 
is violated). FDT is, however, asymptotically recovered 
in the limit A > 1 + r + (which is the limit where c(t) 
follows S(t) closely, see Sec. N} . 

The ratio of cffcctivc-to-real temperature, T (w)/T, 
is shown in Fig. [6] as a function of u for both S and c. We 
see that while the ratio for the channel state approaches 
unity as w — > oo, the corresponding curve for calcium 
has a different asymptotic limit. The system as a whole 
is therefore out of equilibrium already if no feedback is 
present. With feedback, at least the S variable becomes 
equilibrated in the limit of high frequencies. 



1.002 



1.0015 



alto 



1.001 



1 .0005 




(60) 



100 



FIG. 6: The ratio T cB /T (defined in Eq.[60l with the Fourier 
transforms of Eqs. 1471 1481 1541 and I55[l of effective and 'real' 
temperature for the channel state and calcium (inset) is plot- 
ted as a function of the frequency to with and without feedback 
for r+ — 6 and A = 5. Note that the asymptotic limit in the 
case of calcium is different for the two cases and differs from 
1. 



VII. NUMERICAL RESULTS AND DISCUSSION 

We carried out numerical simulations in order to verify 
our analytical results in the weak feedback limit, as well 
as to study the behavior of the system in the intermediate 
(a ~ 1) and strong (a ^> 1) feedback cases. The dimen- 
sionless version of the system described by Eqs. C MS] was 
simulated using a fixed discrete time step At [46| [47| . 

In general, the feedback effects were found to increase 
in significance with larger A. However, very high val- 
ues of A also tend to make the kinetics of c too closely 
tied to that of S (see Fig. @|, and for this reason, we 
performed most of our simulations with the intermedi- 
ate value A = 5 (unless otherwise indicated). The weak- 
feedback (a = 0.1) case was used for detailed comparison 
with the analytical results. We observed from the data 
that the first order perturbation theory works well up 
to a w 0.2. Since the difference between the curves for 
a = and a = 0.1 is small, in the following, the a = 0.1 
curve is not shown in the main figures. Instead, the com- 
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parison between numerically obtained data points (sym- 
bols) and the analytical expressions (lines) is shown in 
the insets of the figures as differences between the cases 
with a = 0.1 and a = 0. The lines in the main figures 
connect numerically obtained data points [48| . 




0.01 0.1 1 10 100 1000 10000 

r+ 

FIG. 7: The numerically obtained open channel fraction (S) 
in the steady state plotted as a function of the opening rate 
r+. Inset: Difference between the cases a — 0.1 and a = 0. 
The line represents the term linear in a from the analytical 
expression in Eq. [40] For the high r+ values as well as for 
a = 100, a smaller time step of At = 10 -4 was used. 

In Fig. the steady state open fraction (S) (which 
is the same as (c)) is plotted as a function of the open- 
ing rate r+, for various values of a. This is analogous 
to a dose-response curve, if the open channel fraction 
(or, equivalently, the mean Ca 2+ ) is interpreted as the 
response of the system to an external stimulus that mod- 
ifies the opening rate of the channels. We observe that 
feedback shifts the response towards higher r + values by 
closing the channels more often. However, the 'dynamic 
range' of sensitivity (e.g. the range of r + covered between 
5% and 95% of the response) is found to be increased with 
feedback: the curves become less steep for large a. 

In Fig. Owe show the variances (5S) 2 and (5c) 2 of S 
and c, respectively, in steady state as a function of r + 
for various feedback strengths. The fluctuations follow 
a bell-shaped curve as a function of r+, but the maxi- 
mum shifts towards larger r + with increasing feedback. 
As we remarked in Sec. El the maximum of the mean 
squared fluctuations (SS) 2 is always | irrespective of the 
feedback strength (and occurs when r + is comparable to 
the average closing rate l + a(c)). For the c fluctuations, 
however, we find that the fluctuations are generally sup- 
pressed by feedback (except at very high r+). The fluc- 
tuations again follow an approximate bell-shaped curve, 
but the maximum sharply comes down with increasing 
a and its position shifts to larger r + . Note also that for 
a > 0, the position of the peak for c-fluctuations always 
trails the corresponding point for channel state. 

To investigate how fluctuations affect the precision of 
signal transduction, we looked at the ratio of the stan- 
dard deviation to the mean (i.e. the noise-to-signal ratio 
in the steady state or the coefficient of variation) for c(t) . 
In Fig. [5] we show 5c/ (c) plotted against r+. For A = 5, 




B r+ 



FIG. 8: Variance in S (A) and c (B) in the steady state 
plotted as functions of r+. Inset: Difference between the 
cases a = 0.1 and a = 0. The line represents the term linear 
in a from the analytical expressions in Eqs. 1511 For the high 
r + values as well as for the case a = 100 a smaller time step 
of At = 10 -4 was used. 



we find that this ratio is a monotonically decreasing func- 
tion of r + , independently of feedback. More surprisingly, 
the ratio itself is enhanced by increasing feedback. Note 
that for r + < 0.1, this ratio is of the order of 10 and 
is reduced to "acceptable" levels (5c/ (c) ~ 1) only for 
high r+ ~ 1 — 10 depending on the feedback strength. 
However, the reverse effect is observed when the pump- 
ing rate is reduced to A = 0.05. In this case, for small r + , 
the noise to signal ratio is reduced by feedback (Fig. [T0|). 

We now turn to time-dependent quantities. Figs. fTTI A 
and B show the auto-correlation functions for the chan- 
nel state and calcium, respectively, in the case of strong 
signal (r+ =6). In general, one notes that the effect of 
increasing feedback is to reduce the characteristic time 
scale of decay of the correlations. When the feedback pa- 
rameter is very large, the channel auto-correlation func- 
tion briefly becomes negative before vanishing at longer 
times. This anti-correlation in the channel state may 
be interpreted physically as the rapid closing of an open 
channel by the Ca 2+ which enters through it. In other 
words, in the presence of feedback, a channel which is 
open at any point of time has a reduced open proba- 
bility at subsequent times (compared to the uncondi- 
tioned mean open probability). By contrast, the auto- 
correlation function for c(t) for the same parameters de- 
cays monotonically with time. In the cases of weak signal 
(r+ = 0.5) or high pump rate (A = 50), both the chan- 
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10 



100 



FIG. 9: The coefficient of variation plotted as a function of 
r+ for A = 5. Inset: Difference between the cases a — 0.1 
and a = 0. The line represents the term linear in a from the 
analytical expression in Eq. [52] For the high r+ values as well 
as for the case a = 100 a smaller time step of At = 10 -4 was 
used. 




100 



FIG. 10: The coefficient of variation as in Fig. [9] but with 
A = 0.05 (only the curves for a = and 100 are shown for 
clarity). Note that, in contrast to the previous case, the ra- 
tio is reduced with increasing feedback for small r+, which 
was also predicted from the first order perturbative correc- 
tion term in Eq. 1521 



nel state and calcium auto-correlation functions decay 
monotonically (data not shown). 

In Fig. Q21 we plot the power spectral density for cal- 
cium and the channel state calculated from the auto- 
correlation function using Eq. [9] in the presence and 
absence of feedback. We observe that irrespective of 



feedback, P c {u>) 



and Ps{oj) 



as u 



It may also be noted that Ps(u>) is non-monotonic in 
this plot, and the observed peak is consistent with the 
dip in the corresponding auto-correlation function in S 
(Fig. [TTA). Notice also that a similar peak is absent for 
P c , which is also consistent with the fact that C c (t) is 
monotonic (Fig. [TIB). In a different context, the reduc- 
tion in noise strength and the shifting of the peak of the 
power-spectrum from low to high frequencies with in- 
creasing feedback has also been observed in simulations 
of a model of an autoregulated gene circuit [29[ with nega- 
tive feedback. An instructive discussion of the the general 
conditions under which feedback can lead to a reduced 
noise intensity and increased noise bandwidth is given in 




FIG. 11: Auto-correlation function for S(t) (A) and c(t) (B) 
in the steady state for A = 5 and r+ = 6. Inset: Difference 
between the cases a = 0.1 and a — 0. The line represents 
the term linear in a from the analytical expressions for Cs{t) 
(Eq.rJH) and C c (t) (Eq.rjSJ. 
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FIG. 12: Power spectral densities for the channel state (main 
figure) and calcium concentration (inset) are plotted against 
the angular frequency ui with and without feedback on a log- 
arithmic scale. The parameter values used are r+ = 6 and 
A = 5. The straight lines are guides to the eye. To get results 
also for high frequencies, a small time step of At = 10~ 5 was 
used. 

Fig. [13] shows the time evolution of the open- 
probability of the channel (Fig. [T31A) and the mean cal- 
cium concentration (Fig. [13)3) for different values of the 
feedback parameter a when the channel was initially 
closed and c(t = 0) = 0. Note that the relaxation to 
the steady state is monotonic for small a, but for higher 
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a, the response in the channel (Fig.[13lM exhibits a max- 
imum, and then settles into the steady state, while the 
response of c does not show this behavior (Fig. [13)3). 




B t 

FIG. 13: Time evolution of the open-probability (A) and the 
mean calcium concentration (B) when the channel starts in 
the closed state (S(t = 0) = c(t — 0) — 0) for different values 
of the feedback parameter a. The opening rate is r+ — 6 
and A = 5. Inset: Difference between the cases a = 0.1 and 
a = 0. The solid line represents the term linear in a from the 
analytical expressions Eq. EH] (0, 0, t)) and Eq. [S] 

In Fig. 1141 we show the linear response function xs(T) 
for different values of a. In order to compute the re- 
sponse function numerically, the system was evolved with 
r + = 0.5 until it reached the steady state. Then, r + was 
increased to 0.6 and the response function was calcu- 
lated using Eq. 1101 By definition, the response function 
shows how a sharply peaked input stimulus (<fi(t) ~ 5(t)) 
is 'transmitted' across the channel. We observe that in- 
creasing the feedback reduces the time constant of decay, 
making the output signal sharper and more similar to 
the input, when A is relatively high. This is in agree- 
ment with our assertion that high values of A improve 
the short-time scale response characteristics of the sys- 
tem. We also confirmed this by studying the effect of A on 
the response time when feedback is strong: in this case, 
reducing A was found to increase the response time ap- 
preciably (for example, when a = 10, a 50-fold reduction 
in A was observed to almost double the response time). 

We used the channel linear response function derived 
in Sec. [VI] to compute the response of the system to 
a periodic stimulus of the form <p(t) — a smut for 14 
different values of to (Fig. [15]) . An explicit calculation 
shows that the linear channel response has the form 




0.5 1 1.5 2 2.5 3 

t 



FIG. 14: The linear response function of the system, com- 
puted numerically using a step stimulus (see text) and a nu- 
merical differentiation for various values of feedback (A = 5). 
The points shown are the averages over 10 9 independent runs 
(symbols are used to indicate the larger error in this calcu- 
lation - the lines in the main figure are guides to the eye). 
Inset: The difference between a — 0.1 and a — is com- 
pared with the analytical expression in Eq. [54] (line) . 

(S(t)) - (S) = Asin(ut + 9), where both the ampli- 
tude A and the phase lag 9 are frequency-dependent: 
A(u) = o/(l + r + )[(l + r+f + w 2 ]" 1 / 2 + 0(a) and 
9(uj) = —2 + arctan[(l + r + )/uj] + 0(a), omitting the 
0(a) terms for the sake of brevity. The complete expres- 
sions (including the first order corrections (to be found 
at [30j)) for A and 9 are checked against numerical simu- 
lations in the insets of Figs. [15] A and B. The amplitude is 
always a decreasing function of the frequency, however, 
the feedback has the effect of increasing the amplitude 
of the response relative to the no-feedback case at suf- 
ficiently high frequencies. At lower frequencies, the op- 
posite effect is observed: introducing feedback tends to 
reduce the amplitude of the output, however, the decay 
as a function of the external frequency is more gradual. 
In simple terms, the over-all effect of the feedback here 
is to widen the range of frequencies effectively sensed by 
the system. The increase in response at large frequencies 
is consistent with the reduction in the time constant of 
the linear response function which we noted earlier. At 
small frequencies, on the other hand, the overall adapta- 
tion produced by the negative feedback leads to a drop 
in the response. 

The phase shift of the response decreases monotoni- 
cally from to — ? as ui is increased, although the change 
is more gradual when feedback is present. Note that for 
large a, there is a significant linear regime in the 9 — to 
curve, which means that in this regime, the time delay 
of the response is effectively independent of the stimulus 
frequency. 

To summarize, numerical simulations provide excel- 
lent support for the analytical predictions from the path- 
integral theory in the weak feedback limit. In addi- 
tion, we observed qualitatively new features in the auto- 
correlation function and in the relaxation of the mean 
open probability to the steady state when feedback is 
strong. In particular, Figs. [51 [TU1[TT1 and [T5l are the prin- 



13 




-0.5 1 1 1 1 1 1 1 1 1 1 1 

2 4 6 8 10 12 14 16 18 20 



FIG. 15: Amplitude (A) and phase (B) of the response (S(i)} 
to a sinusoidal signal (j>(t) = asm(u)t) with amplitude a = 0.1 
and baseline r+ = 1, plotted as a function of stimulus angular 
frequency ui (A = 5). The lines in the main figures are guides 
to the eye. Inset: Difference between the cases a = 0.1 and 
a — 0. The line represents the 0(a) correction computed 
analytically (not shown in text). Note the regime of linear 
decay of versus u in B for small w, which signifies a constant 
time delay in the response. 



cipal results of this paper, and will be discussed more in 
the next section in the context of experimental results. 



VIII. CONCLUSIONS 

In this paper, we studied the stochastic kinetics of a 
simple auto-regulatory signaling module with negative 
feedback using path-integral techniques and numerical 
simulations. Within our formalism, all the statistical 
averages and correlation functions that characterize the 
long-time kinetics of the system are formally expressed 
as a power series in the feedback parameter. Explicit 
expressions were obtained for mean values and correla- 
tion functions to first order in the feedback parameter. 
These analytical results were compared to numerical sim- 
ulations. The simulations were also done in the strong 
feedback regime beyond the applicability of perturbation 
theory. In particular, we investigated how the system 
responds to temporal variations in an external stimulus, 
and how the feedback regulates this response. 

Our principal conclusions from this study may be out- 
lined as follows. We find that, in the context of the ion- 
channel model, the rate of removal of Ca 2+ ions, A, is the 



single most important control parameter of the model. 
When A is large (compared to the intrinsic closing rate 
for the channel), the short-time scale dynamic response 
characteristics of the system are improved by negative 
feedback. In this case, feedback reduces the time scale 
of decay of the auto-correlation and response functions, 
and thus extends the response of the system to higher 
frequencies of input stimulus. 

For all parameter values, the negative feedback shifts 
the range of sensitivity of the system (i.e. the range of 
stimulus magnitudes for which the response is apprecia- 
bly non-zero and not saturated) towards stronger stimuli. 
At the same time, the range of sensitivity gets broader 
(i.e., the "dynamic range" is increased) with increasing 
feedback. While the steady state response to a fixed 
stimulus is always reduced by the negative feedback, a re- 
markable effect was seen in the coefficient of variation for 
the calcium concentration. Increasing the feedback was 
found to increase this quantity in the large A regime. For 
small A, however, the opposite was seen: the coefficient of 
variation was reduced with feedback, for sufficiently weak 
stimuli. In this regime, the negative feedback therefore 
improves the precision with which one can discriminate 
among stimuli of different magnitudes based on the sys- 
tem output. Having a small A, however, was generally 
found to adversely affect the temporal response charac- 
teristics. 

Starting from the experimental literature, we have ob- 
tained estimates for the dimensionless parameters in our 
model for the specific case of calcium signaling in olfac- 
tory sensory neurons (the corresponding module is shown 
in the lower part of Fig. [2]). The rate of opening of 
the cyclic-nucleotide-gated channel varies in the range 
r + i=s 10 -5 - 20 depending on the strength of the odor- 
ant stimulus. The parameter governing the calcium dy- 
namics is A « 20 and for the feedback parameter we find 
a 1 - 10 (4{|. Note, in particular, that A is found 
to have a relatively large value. In the view of our re- 
sults, this suggests that the calcium-mediated feedback 
in olfactory transduction evolved so as to improve the 
temporal response characteristics of the system, rather 
than to improve the discriminability of weak stimuli. 

Our findings are consistent with the results of ex- 
perimental studies in a variety of sensory systems with 
calcium-mediated negative feedback. The feedback- 
induced shift of the range of sensitivity to higher stimu- 
lus magnitudes has been demonstrated in olfactory sen- 
sory neurons and in cone photoreceptor cells (reviewed 
in Ref. jU S2|). The temporal characteristic of response 
was found to be altered in several studies in which the 
strength of feedback was experimentally manipulated. In 
Ref. [321 ] . the time course of response to a pulse stimu- 
lus was found to be slowed down when the frog olfactory 
cilium was put in lowered external calcium concentra- 
tions (in such conditions, the ion flux J, and therefore 
the dimensionless coupling a, is decreased). A more pro- 
nounced effect of this type was observed in similar exper- 
iments for Drosophila photoreceptor cells in Ref. [31| • A 
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marked slowing down of the response to a step stimulus 
was observed for the newt olfactory cilium in Ref. [43[ 
when the calcium flux J was decreased by changing the 
holding potential. Fluctuations in the steady state have 
received less attention in the biological literature. The 
power spectrum of odorant-induced current fluctuations 
across the membrane of rat olfactory sensory cilia was 
studied experimentally by Lowe and Gold [33j | . They ob- 
served that the tail of the power-spectrum decays with 
an effective exponent in the range 2.3 - 2.5, which is sug- 
gestive of a crossover between uj~ 2 and uj~ a decay. In 
this case, however, the dynamics of the cAMP module 
(shown in the upper part of Fig. [2]) should be taken into 
account in the theoretical treatment, as fluctuations in 
cAMP give an important contribution to the variability 
of the response [33| ■ 

We are currently in the process of extending this study 



to include spatial effects. In particular, in olfactory sen- 
sory neurons, the Ca 2+ ion channels are spatially dis- 
tributed along long and thin cellular compartments called 
cilia. The feedback studied in the present manuscript 
coupled with calcium diffusion inside the cilium can give 
rise to non-trivial spatial correlations between the chan- 
nels, and it will be interesting to study how the spatial 
coupling affects the overall kinetics of the system. 
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[45] 
[46] 



[47] 



[48] 



[49] 



/ T>{T} n Pu[t n -2, t n -V, {r} n ; c„_ 2 ; 2m„]. In the presence 
of a feedback, the weights Vij[to,t;{Tk}k=i',co; N] are 
functions of the concentration cq at time to. Since for 
any realization S(t), the function c(t) follows according 
to Eq. [5j the values Ck = c(tk) are nontrivial. As a con- 
sequence, correlation functions of order larger than two 
cannot be written as simple products of propagators. 
Putting in the coefficients from Appendix iDl shows that 
the terms ~ lo~ 2 cancel also for a ^ 0. 
The time step was chosen to be much smaller than all the 
relevant timescales: At <C mini—, tt— , tI. The kinetics 
of the (dimensionless) calcium concentration is computed 
using a simple Euler forward step algorithm. (Unless oth- 
erwise indicated, we used At — 10 -3 and the curves pre- 
sented are averages over 10 6 independent runs.) Feedback 
strengths a = 0, 0.1, 1, 10 and 100 were used, and r+ val- 
ues were chosen from the range [0.01, 10000]. 
In principle, a Gillespie-type algorithm [23] could be 
used. Its implementation would, however, require numer- 
ically inverting the dwelling time distribution Pi(t) = 
exp{-i [a(c(r 3 -) - 1)](1 - e- xt ) - (1 +a)t} for the open 
state for each value of c (see also [Sj). It is not clear, if 
this approach would decrease computer time compared 
to the fixed time step algorithm. 

The data points in the main figures (many hundreds in 
the plots over time, approximately 200 in the plots over 
r+) are equally spaced on the linear or logarithmic scale. 
Through variation of At and the number of independent 
runs used to calculate the averages, it was checked, that 
the systematic error due to time-discretization and re- 
stricted sampling size is in the range of line thickness. 
In the insets as well as in the main figures of Figs. [14] 
and 1151 this error is smaller than or comparable to the 
symbol size. 

To obtain these estimates, we view the olfactory cilium 
as consisting of a series of cylindrical segments, each 
containing one cyclic-nucleotide-gated (CNG) channel. 
On the time scales of interest, we neglect calcium dif- 
fusion between neighboring segments. We estimate the 
parameter values in each segment as follows: J « 1.3 ■ 
10" 19 mol ■ s _1 V » 6 • 10" 24 m 3 (calculated from 

the density of CNG channels and the geometry of the 
cilium as cited in [23], assuming a homogeneous distribu- 
tion of channels); A ~ 10 4 s _1 (calculated from a Taylor 
expansion of the calcium-extrusion term used in the de- 
terministic model of Ref. [22J); Rl_ » 500s _1 (from single 
channel measurements with weak stimulus and low cal- 



cium (e.g. [35j, [3a])); R+ ps 5 • 10 



(based on 



an open probability of 10 if no stimulus is present to 
0.95 for full activation gj}); a » 3 • 10 8 - 2 • 10 9 s" 1 M" 1 
(based on data in the figures of Refs. [3^] (lower estimate) 
and [13] (upper estimate)). 



APPENDIX A: CALCULATION OF INTEGRALS 
7 , Ii AND I 2 

Generalized convolution theorem: 

The following general result is very useful in perform- 
ing many calculations using the path-integral technique. 

if g(t) = Si Hh) Jo"* 1 h(t2) ■ ■ ■ /*"•-*»-* f n (t n ), then 



its Laplace transform g(s) = / °° g(t)e st dt is given by 

n 

g(s) = -l[f~M- (Ai) 



This result is proved in a straightforward way by re- 
peated application of the standard convolution theorem 
of Laplace transforms. 
Calculation of Iq: 

The integral Iq is relatively easy to compute. The 
Laplace transform Iq(s) can be written easily from the 
definition of the integral in Eq. [M] and using the above 
theorem: 



I (s;m) = s- (m+1 )(s + l-r + )- 



(A2) 



Calculation of 1\\ 

In order to calculate the integrals I\ and I2 , it is neces- 
sary to express c(t) in terms of the time interval variables 
using Eq. [5j The explicit relations are, 



c(t) =1 - e - x( - T - Tl) for ti<t<t 2 , 

c(t) =1 + e" Ar (e Ar2 - e Ari + • • • + e^- 1 - e AT >) 

for Tj < t < Tj+i with odd j > 3. (A3) 

From that it follows 

i— 1 

z(n) = [e- x(T ^+ l) - e-^--^)) for even i > 2. 



i=i 



(A4) 



Using Eq. IA4I in Eq. [MJ we write I\ — g\ — #27 with 

2m i — 1 



91 



(0,i;m)=]T' 5^" y VT 



«=2 j=l 



exp 



(l-r + )53*i+A(7t-T i+1 ) 



k=l 



(A5) 



and 



2m i—1 



.92 



(0,«;m)=53' J2" I VT 

A — O A — 1 J 



1=2 3=1 



exp 



(l-r+^t'k + Kn-Tj) 



k=l 



(A6) 



Note that when j < i-3,Ti—Tj + i = Y,%+s\/2 (*fc + ife ) 
and is zero when j = Let us therefore split the 

second sum, and treat these cases separately. In order to 
compute the Laplace transforms, we use the theorem in 
Eq. IAU It is now straightforward to write the Laplace 
transform of gi in Eq. IA5l using this theorem. The result 
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where 



«h(s) = ms-( m+1 \s + 1 - r+y 



2m i — 3 

/ . ^ — ir 



-0'+l)/2 



i=4 j=l 



(A7) 



where we have defined y = s(s + A) _1 (s + 1 — r+)(s + 1 — 
r + +A) _1 < 1. The double geometric sum in Eq. lATl turiis 
out to be j/) -1 [m + (1 — y)~ 1 (y m — 1)] , which gives 

- ,s- (m+1) ( S + 1 - r+)- m -^x 

1-2/ 



m + -^(y m -l) . (A8) 
1-2/ 



Calculation of c/2 also proceeds along similar lines. Af- 
ter some straightforward algebra, we find that 

aw =.-<-+■■(.+ i-r t) ~ ,'^; r ;, » 

2m i/2-1 \ 

<" + £V /2 £ w - ' • ( A9 ) 

j=2 ;=i y 

The double sum is easily shown to be equal to y(l — 
y) _1 [rn + (1 — y)^ 1 {y m — 1)], and this gives 

h{s) = S - (m+1) (, + 1 - r + )-" ~j~~~~~~T x 

s + 1 - r+ + A 

' ^m+-^-( tf m -l)V (A10) 



1 - 2/ V 1-2/ 
After putting together Eq. IA8I and Eq. IA101 we arrive at 

s + X 



Ii(s;m) — s~( m+lS) (s + 1 - r + y 



2s + 1 - r + + A 



y 



m+-^(y m -l) . (All) 
1-2/ 



Calculation of li : 

For ^2 , the following integral is needed (obtained using 
Eqs.[53J): 



1 



C (r)dr = ^ +1)/2 + -(e-^-e- A ^) 



j'-i 



5] - £ ^< (A12) 



fc=2 



;=i 



for odd j > 1. 

The integral I2 may be expressed in the form 

l 2 (0,t;m) =h + h 1 + -(h 2 + h 3 -hi- h 5 ), (A13) 

A 



(r+ 



/to = / PTexp 

2m-l 

fti = J^" / X>T exp 

i=3 J 

2m-l j-l 

^=£" £'/ PTex p 

j=3 fc=2 ^ 
2m -1 j 

ft 3=E" E" / PTex p 

3=3 1=1 ^ 
2m -1 j 

^ 4 =E" E" / PTex p 

3=3 1=1 ^ 
2m- 1 j-l 

/i s=E" E' / PTex p 

V ^ f o « 



fc=l 



c{t)cIt 



t' 



0'+l)/2 



(r + -l)£4 
k=l 

( r +- 1 )E^ + A (^- Tfc 

m 

(r + -l)^/' fc + A(r J+1 -T / ) 
fc=i 

m 

( r+ -l)][Y fc + A(7- i -7i) 



fc=i 



j=3 fe=2 



( r + - i)E^ + a ( t j+! ~ Tk 

1=1 

(A14) 



Obviously, ft represents j = 1 in the sum in Eq. 1341 and 
hence includes all possible contributions from m = 1. 
The rest of the terms therefore always have m > 2. The 
Laplace transforms for ho and hi are straightforward, 
and one may easily verify that 



ho(s) =s-( m+1 \s + 1 - r+y m x 
1 1 



s + 1 — r+ s + 1 — ?*+ + A 
m — 1 



ft! (a) = s -(" l+1 )( s + l-r + )- 



s + 1 — rj 



m > 2 
(A15) 



so that 



ho + k =s-( m+1 \s + l-r+y 



s + 1 — r + s + 1 — r_j_ + A 



(A16) 



for m > 2. We omit further details of computation of 
^2,^3,^4 and which involves straightforward, but 
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somewhat tedious algebra. The results are (for m > 2) 
y m - 1 \ s s + l-r + +X 



1-y J\2s + l-r + +\ 



h(s) =s- (ro+1) (s + l-r+) 
, Z/(y m -l) 



(l-*))(4 + S + 1 ^ + 



1-y v °'J\\ ) 2s+l-r+ + A 
fc 4 (s) = s _(m+1) (s + 1 - r+)- m x 

y(y m - 1) ,„ A /s , „\ s + l-r+ + A 



(1 -»>) G + 1 )STT 



^s(s) = s~ (m+1) (« + 1 - r+)- m x 

y m — 1 \ s s + 1 — r_| 



1-y / A 2s + 1 - r + + A 



- r+ + A 



(A17) 



After putting together all the terms, we show that 
I 2 (s;l) ^(s+l-r+y 1 

(A18) 



1 



s + 1 — r+ s + 1 — r+ + A 

and for m > 2, 

J 2 ( S ; m ) = S -(" l+1 )(s + 1 - r+)- m x 

m(s + A) 



+ 



(s + 1 - r+)(2s + 1 - r+ + A) 
y m - 1 s 
1-y (s + 1 - r+ + A)(2s + 1 - r+ + A) 



(A19) 



APPENDIX B: RELATION BETWEEN Gn AND 
Goo FOR a 7^ 

In order to proof Eq. [37j to be valid up to 0(a), 
we first calculate the functions fi(t) and /2(f)- As in 
Eq. EU in Eq. EH c'/i(i) and c'/ 2 (t) are the 0(a) terms 
in Goi(0, c';i) and Gu(0,c';t), respectively, when the 
channel does not allow entry of ions. The solution is ob- 
tained by solving the corresponding rate equations with 
the time-dependent rate r-(t) = l + ac'e~ xt . The results 
are 



/i(t) =-/(*). 



-(l+r+)t _ 



A V 1 + r+ - A 
(l + r+)(l + r+-A)' 



1 + r+ 



,-At 



-(l+r++A)t 



(Bl) 



full propagators II with < i' < t: 

G i(0,0;i) = V / tic'/ dcn 0;7 -(0,0;t / ,c')II J -i(t',c';t ) c) 
„• Jo Jo 



dc'n O o(0,0;i',c')Goi(^,c';t) 



/ dc'Il 01 (0,0;l/,c')Gn(t',c';t). (B2) 
Jo 



After substituting Eq. [39] into Eq. IB21 we find that 

G01 (0, c = 0; t) = G oo (0, c = 0; t')G 01 (t\ c' = 0; t) 
+ G i(0,c = 0;tOG ! ii(f , ,c' = 0;t) 

+ a (/!(« - £')<</ (t')>o + / 2 (t - 0<C(0>i) + 0(« 2 ), 

(B3) 

where (c'(f'))o and (c'(i'))i were defined in Eq. [13] Let 
us now take the limit of f' — ► as in the case without 
feedback, and express the Green's function Gn in terms 
of Goo- The result, to 0(a) is 

G 1 i(0,c o = 0;t) = l-Goo(0,0;f) 
1 



[d t G oo (0,0;t) 



a t G oo (0,0;t)| t=o 

-a(/i(t)^(c'(t'))o|t'=o + /2(t)9 4 '(c / (t / ))i| t '=o)]. 

(B4) 

It is easy to see now that both the time derivatives 
of (c'(£')) /i in the previous equation vanish, simply be- 
cause they are defined in Eq. [13] with an initial condition 
c(0) = by construction. This is easily seen by first solv- 
ing Eq. [5] to find (c(t)) to order a = 0, and the result 
is 



(c(t)) 



1 



(l + r+)\ l + r+-A 
(1 + r+)e- xt - Ae- (1+r +H J + 0(a), (B5) 



Instead of Eq. [35] we now write a relation between the 



whose time derivative vanishes at t — 0. But, since 
(c(t)) = (c(t)) + (c(i))i and since ft{c(i)Wi |t=o > 0, 
it follows that <9 t (c(i)) /i |t=o = in Eq. IB4I and this 
leads to Eq. [571 



APPENDIX C: CALCULATION OF (c) AND (c)i 

In this appendix, we compute (c(£))o, whose long-time 
limit gives (c)o- Note that, in accordance with Eq. IB31 
this quantity needs to be computed only to the zeroth 
order in a. This can be done using the probability func- 
tional Voo from Eq. [23] We therefore define 

OO r. 

(c(t)) =53/ PT7, oo[0, i; {*<}, {tj}; 2m]c(i), (CI) 
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//2m— 1 



where c(t) — e xt ^ 



from Eq. [5] 

Further calculations are done with the aid of Laplace 
transforms and using the generalized convolution theo- 
rem (Eq. lAlj) . We omit further details, and simply give 
the final result for the Laplace transform (c(s))q: 



(c(s))o 



r + X 



s(s + A)(s + 1 + r+)(s + 1 + r+ + A) ' 



(C2) 



The steady state value is obtained after inversion s — ► t 
and taking the limit t — > oo and is easily found to be 



Co = 



r+ 1 
1 + r+ 1 + r+ + A ' 



from which we also find 



c 1= c - c 



1 + r+ 1 + r+ + A ' 



(C3) 



(C4) 



B 9 



r+A 



C 2 = 



(1 + r+) 3 (l + r+ - A) 3 (l + r+ + A) 2 
[(A - l)rj + (A - 1)(2 - A)r 3 - A(A 2 - 4A + l)r£ 
+ (A 4 -A 3 + 7A 2 -7A + 2)r + 
-(A-1)(A 3 - A 2 -3A + 1)] 

r J± > 

(1 + r+) 3 (l + r+ - A) 3 (l + r+ + A) 2 (l + r+ + 2A) 

[r£ + (1 + X)r% - (3A 2 - 5A + 6)r* 



+ (-A 3 +4A 2 + 6A-14)4 

+ (2A 4 + A 3 + 18A 2 - 2A - ll)r 2 



D 2 

Fn = - 



+ (A 4 


+ 3A 3 + 12A 2 - 7A - 3)r+ 


+ A(A 3 




r+A(r+ - 1) 




(1- 


f r+) 3 (l + r+ + A)(l+r + 


+ 2A) 




r + A 2 (A- 1) 




(1- 


f r+) 3 (l + r+ - A) 3 (l + r_ 


f + A) 2 


[r% + 


3r 2 + + (3 - A 2 )r + - (A 2 - 


1)] 




r 2 A 2 




(1- 


f r+) 2 (l + r+ - A) 2 (l + r_ 


f+A) 



APPENDIX D: LIST OF COEFFICIENTS IN THE 
FIRST ORDER TERMS 



1. Correlation functions 



2. Response functions 



Bi = 



Ci - 
Di = 
E x = 



(1 + r+) 3 (l + r+ - A) 2 A(1 + r+ + A) 
[-(A - l)r% - (A - 1)(2 - A)r 3 + A(A 2 - 2A - l)r\ 
+ (-A 4 + A 3 - 3A 2 + 3A - 2)r+ 
+ (A-1) 2 (A + 1)(A-1)] 
2r 2 A 



(l + r+) 2 (l + r+ -A) 2 (l + r+ + A) 
r + (ri-l) 



(l + r+) 3 A(l + r+ + A) 

r+(A-l) 
(l+r + ) 2 (l + r+-A) 



S 3 



1 



-(A-l)r 



C 3 



r> 3 = 



#3 = 



A(l + r + ) 2 (l + r + — A) 2 
+ (2A 2 - 4A + l)r\ + (-A 3 + 2A 2 - A - 1)? 

-(A -I) 2 ] 

r+X 

(l + r+) 2 (l + r+ - A) 2 
rl-X-1 



(l + r+) 2 A(l + r+ + A) 
A- 1 

(l + r + )(l + r+-A) 
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B4 = - (l + r + )*(l + r + -A)3 [1 + ^-3) 



C 4 = 



£>4 

F 4 



+r+(l + A 2 (A - 1)) - r\ (1 + 2A(A - 2)) + r^(A - 1) ] 

A 

(l + r+) 3 (l + r+-A) 3 (l + r + +A) X 
[A 3 (l + r++r 2 ) 

-A 2 r+(1 + r+) 2 - A(l + r+) 2 (l + r+(r+ - 5)) 
+r + (r + -l)(l + r+) 3 ] 

- A - 1 
~ (l + r+) 3 (l + r+ + A) 
A(A-l) 
(l + r + )(l + r+-A)2 

r+A 2 



(l + r+)2(l+r + -A)2 



